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The  shell  theory  used  to  model  each  conical  frustum  is  a 
third  order  Flugge-Lure-Byrne  theory  in  which  a complete  inter- 
action between  shell  membrane  and  bending  effects  occurs. 

Energy  forms  are  used  to  generate  stiffness,  mass  and  x\ 

stability  matrices  required  for  the  finite  element.  The  stability 
functional  employed  is  Bolotin's  initial  stress  functional  to 
which  is  appended  an  external  potential  for  normal  pressure 
loadings  as  defined  by  Cohen. 

All  variables  used  to  describe  the  system  are  represented  by 
means  of  Fourier  series  in  the  circumferential  angle  0.  Except 
for  the  stability  functional,  a complete  decoupling  exists  so 
that  only  terms  associated  with  a given  circumferential ‘ mode 
number  n interact.  For  the  stability  functional  only  the  initial 
stresses  associated  with  the  symmetric  pressure  field  are  con- 
sidered in  order  to  complete  the  Fourier  decomposition,  but  both 
initial  membrane  and  bending  stresses  are  retained. 

The  kinematical  behavior  of  the  displacement  field  of  the 
conical  frustum  segment  associated  with  each  mode  number  n is 
defined  by  two  different  models.  For  the  first  (simpler)  model, 
the  normal  (outward)  displacements  are  third  order  polynomial 
functions  of  the  position  along  the  generator,  while  the  in-surface 
displacement  fields  are  first  order.  For  the  second  (refined) 
model,  the  normal  displacements  are  fifth  order  polynomials  and 
the  in-surface  displacements  are  third  order.  Thus  the  simpler 
model  of  the  conical  frustum  has  two  nodal  planes  (perpendicular 
to  the  axis  of  the  frustum)  and  four  degrees  of  freedom  at  each 
node,  totalling  eight  degrees  of  freedom  for  the  element.  The 
refined  element  has  three  nodal  planes  again  with  four  degrees  of 
freedom  at  each,  so  that  the  finite  element  model  has  twelve 
degrees  of  freedom.  All  integrations  required  to  generate  the 
stiffness,  mass  and  stability  matrices,  are  performed  numerically, 
using  Gaussian  quadrature. 

The  finite  element  models  presented  herein  are  to  be  basic 
elements  in  an  anlysis  program  designed  specifically  for  use  in 
a program  for  the  automated  structural  design  of  submerged  stiffened 
pressure  vessels  in  the  form  of  arbitrary  structures  of  revolution . . 
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ABSTRACT 

This  report  presents  the  analytical  development  of  a finite 
element  model  suitable  for  a study  of  the  elastic  vibration  and 
stability  characteristics  of  a structure  composed  of  a finite 
number  of  shells  of  revolution,  each  in  the  shape  of  a conical 
frustum.  Thus,  the  model  can  be  used  to  study  a shell  of  revo- 
lution, with  arbitrary  meridional  cross  section,  to  which  is 
attached  arbitrarily  shaped  circumferential  ring  stiffeners  and 
lateral  bulkheads. 

The  shell  theory  used  to  model  each  conical  frustum  is  a 
third  order  Flugge-Lure-By rne  theory  in  which  a complete  inter- 
action between  shell  membrane  and  bending  effects  occurs. 

Energy  forms  are  used  to  generate  stiffness,  mass  and  stability 
matrices  required  for  the  finite  element.  The  stability  functional 
employed  is  Bolotin's  initial  stress  functional  to  which  is  appended 
an  external  potential  for  normal  pressure  loadings  as  defined  by 
Cohen . 

All  variables  used  to  describe  the  system  are  represented  by 
means  of  Fourier  series  in  the  circumferential  angle  Q.  Except 
for  the  stability  functional,  a complete  decoupling  exists  so 

I 

that  only  terms  associated  with  a given  circumferential  mode 
number  n interact.  For  the  stability  functional  only  the  initial 
stresses  associated  with  the  symmetric  pressure  field  are  considered 
in  order  to  complete  the  Fourier  decomposition,  but  both  initial 
membrane  and  bending  stresses  are  retained. 

The  kinematical  behavior  of  the  displacement  field  of  the 
conical  frustum  segment  associated  with  each  mode  number  n is 
defined  by  two  different  models.  For  the  first  (simpler)  model, 

i 
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the  normal  (outward)  displacements  are  third  order  polynomial 
functions  of  the  position  along  the  generator,  while  the  in-surface 
displacement  fields  are  first  order.  For  the  second  (refined) 
model,  the  normal  displacements  are  fifth  order  polynomials  and 
the  in-surface  displacements  are  third  order.  Thus  the  simpler 
model  of  the  conical  frustum  has  two  nodal  planes  (perpendicular 
to  the  axis  of  the  frustum)  and  four  degrees  of  freedom  at  each 
node,  totalling  eight  degrees  of  freedom  for  the  element.  The 
refined  element  has  three  nodal  planes  again  with  four  degrees  of 
freedom  at  each,  so  that  the  finite  element  model  has  twelve 
degrees  of  freedom.  All  integrations  required  to  generate  the 
stiffness,  mass  and  stability  matrices,  are  performed  numerically, 
using  Gaussian  quadrature. 

The  finite  element  models  presented  herein  are  to  be  basic 
elements  in  an  analysis  program  designed  specifically  for  use  in 
a program  for  the  automated  structural  design  of  submerged  stiffened 
pressure  vessels  in  the  form  of  arbitrary  structures  of  revolution. 


I. 


INTRODUCTION 


This  report  presents  the  analytical  development  of  a 
finite  element  model  suitable  for  the  analysis  of  an  arbitrary 
elastic  shell-type  structure  of  revolution  subjected  to  arbitrary 
applied  dynamical  loadings  and  an  external  pressure.  In  this 
model,  the  structure  is  composed  of  a finite  number  of  shells  of  revo- 
lution, each  in  the  shape  of  a conical  frustum,  so  the  the  model 

can  be  used  to  study  an  arbirary  shell  of  revolution  to  which  is 
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attached  arbitrarily  shaped  circumferential  ring  stiffeners  and 
lateral  bulkheads.  This  model  and  the  computer  program  based  on 
the  model  is  intended  to  be  used  as  an  analysis  subprogram  in  an 
automated  structural  design  algorithm. 

Computer  programs  for  the  analysis  of  complex  assemblies  of 
shells  of  revolution  have  been  in  existence  for  nearly  a decade 
[1.2].  and  in  particular  the  BOSOR  computer  code  [2] , has  received 
wide  spread  acceptance  and  usage  in  the  engineering  community. 

Codes  of  this  type  are  necessarily  developed  for  user  ease  and 
convenience,  with  special  features  designed  (a)  to  ease  preparation 
and  checking  of  computer  models  of  the  structure  and  (b)  to  provide 
comprehensive  data  manipulation  and  display  features  to  assist  the 
user  in  interpreting  results.  Unfortunately,  these  features  make 
the  codes  very  difficult  to  employ  in  the  development  of  an  auto- 
mated structural  synthesis  program.  In  such  a program,  the  analysis 
of  a shell  of  revolution  is  a subprogram  (or  subroutine)  which 
must  be  called  and  exercised  many  times  during  the  course  of  the 
automated  structural  design  procedure.  Thus  the  analysis  subprogram 
must  be  as  efficient  as  possible  and  be  written  to  minimize  computer 
storage  requirements,  with  no  user  convenience  packages. 


The  difficulties  in  developing  such  an  analysis  subprogram 
are  evident,  since  to  date  design  optimization  algorithms  [3,4,5] 
have  been  capable  of  dealing  only  with  simple  structural  config- 
urations such  as  cylindrical  shells  stiffened  by  only  simple 
T-ring  stiffeners.  Analyses  have  been  based  upon  equivalent 
orthotropic  shell  models,  which  are  unable  to  account  for  the 
effects  of  discrete  stiffeners. 

The  finite  element  model  for  the  conical  frustum  presented 
in  this  report  is  based  on  Flugge-Lure-Byrne  third  order  shell 
thepry  [6],  in  which  a coupling  exists  in  the  shell  constitutive 
relations  between  bending  and  membrane  effects.  Stiffness  and 
kinetic  energy  functionals  are  used  to  derive  "consistent" 
stiffness  and  mass  matrices.  Bolotin's  functional  for  initial 
stress  [7]  and  Cohen's  representation  for  the  potential  due  to 
external  pressure  [8]  are  used  to  generate  the  so  called  "geometric 
stiffness"  matrices  required  for  the  study  of  the  effects  of 
external  pressure  on  structural  stability.  All  variables  in  the 
system  are  represented  by  means  of  a Fourier  series  in  the  circum- 
ferential angle  Q . Except  for  the  stability  functional  a complete 
decoupling  exists  between  variables  with  differing  circumferential 
mode  numbers  so  that  only  variables  associated  with  a given  cir- 
cumferential mode  number  n interact.  In  order  to  decompose  the 
stability  functional  the  initial  stress  state  of  the  structure 
is  assumed  to  be  symmetric,  i.e.,  only  the  initial  stresses 
associated  with  n=0  are  considered,  with  both  membrane  and 
bending  stress  resultants  taken  into  account. 

The  kinematical  behavior  of  each  Fourier  component  of  the 
displacement  field  of  the  conical  frustum  is  defined  by  a poly- 
nomial function  of  position  along  the  shell  generator.  Two 
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different  models  are  defined  based  on  the  complexity  of  this 
assumed  displacement  field.  For  the  first,  simpler  model,  the 
normal  (outward)  displacement  is  represented  by  a third  order 
polynomial  and  the  in  surface  displacements  are  of  first  order. 
Therefore,  this  model  is  similar  to  that  in  [l]  although  the 
element  in  [l}  is  based  on  a Donne 1 1-Mush t ar i-V las ov  theory. 

This  model  results  in  a nodal  ring  at  each  end  of  the  conical 
frustum  around  which  three  nodal  displacements  and  one  rotation 
are  defined.  Thus,  this  element  has  eight  nodal  displacement 
degrees-of-freedom. 

The  second,  more  complex,  element  utilizes  a fifth  order 
polynomial  representation  for  the  normal  displacement  and  second 
order  polynomials  for  the  surface  displacements,  resulting  in  a 
model  with  three  nodal  rings  and  twelve  nodal  displacement 
degrees-of-freedom. 

Integration  of  the  energy  forms  required  for  determination 
of  the  coefficients  of  the  stiffness,  mass  and  stability  matrices 
is  accomplished  by  use  of  Gaussian  quadrature  formulas,  using  three 
through  ten  point  interpolation  formulas. 

Section  II  describes  the  general  analytical  development  of 
the  finite  element  model  and  Section  III  gives  the  matrix  forms 
which  describe  the  stiffness,  mass  and  stability  (initial  stress) 
characteristics  for  the  element.  Section  IV  is  a brief  description 
of  the  computer  code  required  to  generate  the  element  properties 
and  of  some  of  the  sample  problems  run  to  date. 
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II.  ANALYTICAL  DEVELOPMENT 


r? 


II. 1 Model  of  Conical  Frustum. 

In  order  to  model  a complex  shell  of  revolution,  it  is 

first  necessary  to  develop  an  analytical  representation  for  a 

frustum  of  a conical  shell,  see  Figure  1.  For  this  frustum 

let  0 denote  the  angle  between  the  (straight)  shell  generator 

and  its  axis  of  symmetry,  s denote  the  distance  along  the  generator, 

(assumed  to  be  the  middle  surface  of  the  shell).  Then  -s  <s<s  , 

o o 

where  2sq  is  the  length  of  the  generator  of  the  shell.  Also,  let 

r.+  r,  r _ — r i s 

0 denote  the  circumferential  angle  and  - H - — • = r 

2 2 so 

be  the  distance  from  the  generator  at  s to  the  axis  of  the  shell, 

where  r and  r are  the  radii  of  the  circles  at  s = - s and  s = s , 
12  o o 

respectively.  For  a conical  shell,  the  principal  radii  of  curvature 

for  the  middle  surface  of  the  shell  are  R = r sec  <p  and  R = oo  in 

1 2 

the  circumferential  and  meridional  direction,  respectively. 

Finally,  let  £ be  the  ordinate  perpendicular  to  the  shell  surface 
as  shown  in  Figure  1 and  let  t denote  time.  Then,  assuming  that 
the  shell  obeys  Love's  postulate  0]  the  complete  behavior  of  the 
shell  is  described  in  terms  of  middle  surface  displacements  u(0,s,t), 
v(#,s,t)  and  w(#,s,t)  acting  in  the  meridional,  circumferential 
and  normal  directions  to  the  shell  surface  (positive  outward), 
respectively.  That  is,  all  points  in  the  shell  move  according  to 
the  relations  [6]  . 


V v 


(i  + 7 cos^)-  7 w,. 


u = u - fw, 
s s 


v w 


(la) 

(lb) 

(lc) 
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where  u^,  u , and  u^.  are  motions  of  a point  at  s , 0 , £ in  the 
circumferential,  meridional  and  normal  directions,  respectively. 

In  order  to  develop  expressions  for  the  internal  energies 
associated  with  initial  and  elastic  effects,  gradients  of  the 
displacement  field  are  required,  written  in  orthogonal  curvi- 
linear conical  shell  coordinates.  The  components  associated 
with  gradient  of  the  ith  displacement  in  the  jth  direction. 


written  y u , are  as  follows  [9]  , written  as  a 3 3 matrix, 

i 3 


>y  u y u y '1 

‘11  12  13 


V u 

2 1 


V u 

2 2 


V U 
2 3 


V u y u V u 
3 1 3 2 3 3 


(u  +u  sin0+u  cos <t> 


0,0  s 


r 


— u cosi 

r\  •«  6 > 


0 , s 


~1 


, s 


0,t 


S,f 


(2) 


where  R = r + £ cos 4>. 

Note  that  the  components  of  the  symmetric  portion  of  this 


matrix  are  the  local  strain  components  e , 

ij 


6 

€ „ 

5 _ 

11 

12 

13 

e 

6 

£ 

21 

2 2 

23 

e 

f 

5 

31 

c3  2 

33 

1 /u  +u  s in0+u  c o s 1^  /u  +AU  -us  s in0\  2.  /u  - ^Ocosft^ 

R \ 6,0  S £ J2[e,s  R S,0  R J 2 l 0,£  R R J 


s , s 


1 /u  +u  \ 

2 \ s>£  £>SJ 


s y mme  trie 


to 


(3) 


The  strains  associated  with  row  3 and  column  3 are  zero  as  a 
consequence  of  Eq.  (1). 


II. 2 Strain  Energy  of  Shell 


The  strain  energy  per  unit  surface  area  U of  a thin 

o 

isotropic  elastic  shell  is  given  by  the  integral 


U = 
o 


2 

2(1-1/)  6 

0s 


(4) 


•U+£  cos<f>)  <f£ 
r 

where  E is  Young's  Modulus,  V is  Poisson's  ratio  anti  h is  the 
shell  thickness. 

Substitution  of  Eq.  (3)  into  (4),  with  the  aid  of  Eq . (1) 

gives,  according  to  the  3rd  order  Flugge-Lure-Byrne  theory  M 

M 

U = 1^  N £ + N £ + ( N , +N  -_0s  ) e 
° 2 6 6 ss  6 s s6  R 6s 

L l 

(5) 

+ M K +M  K + (M  +M  ) K 
6 6 s s 0s  s6  sd 

where 

£ = l(u,  +u  sin0  + wc  o s0) 

6 7 6 


£ =u,  (6) 

s s 

e ~e  = 1l(Au,  +v,  -sin0‘  v) 

03  s0  2 r d s 


are  the  shell  surface  (membrane)  strains,  and 

K = ri  cos<+>  v,  - JL  w,  - sin0  w,  "] 

0 I r r ' 6 r 66  s 

K - - w , 

s s s 

K = k - cos 0 ( v , - sin0  v ) +_1  ( sin <j>  w , - w , ) 
Os  sd  r s r r r ® s 0 


(7) 


are  changes  in  the  surface  curvature  of  the  shell.  The  symbols 
N and  M i , j = 0, s are  the  shell  membrane  stress  resultants 

ij  ij 

and  shell  bending  moments,  respectively,  and  are  related  to  the 
strain  measures,  ( and  k according  to  the  relations 

y y 


N0  = C 6d+  V€ 


+ D / cos# 


e Q c os#  _ Kq 


N = c [ e + V€a  ) + D / cos#  \ K 


No  = C ( 1-  V ) • e. 


+ D(l-y)  cos#  / cos#  eQ  - Kq 
r r os  — 


N a = C (1-  V)  • € + D(l-y)  cos#  . * 


M0  = ° ( Kd  + *"s  * 60 


Ms  - D “s  + "*0 


S X 


M0S  = D<X-^  (KQs  - e0s  ^ 


Ms0  = D (1-  I')  Ksd 


z 3 / 2 

where  C = Eh/(l-y  ) and  D = Eh  //12(l-I'  ) 


Note  that 


N a + — ^ = N 
Os  d s 0 
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II. 3 Kinetic  Energy 

The  kinetic  energy  of  the  shell  per  unit  surface  area  is 
given  by  the  integral 


/h/2  2 2 2 

(u  + u + u )(1  + Ceos  <f>)  d£ 
h / 2 o s C r 


(10) 


In  (10)  the  dot  denotes  differentiation  with  respect  of  time  t. 

Since  u , u and  u are  related  to  the  surface  displacements 
0 s C 

v,u  and  w according  to  Eqs.  (1),  substitution  of  Eqs.  (1)  into 
Eq.  (10)  gives 


T = ph 
o 2 


< 2 2.2  2 2 
1 + h cos  v + u + w 

' *r'  ) 


- h cos0  ./  2vw,  + ruw, 

2 • \ 0 s 

, ' 1 
6r 


+ h 


2 

w», 


12  \ 2 0 


2 

+ w, 


(ID 


II. 4 Potential  for  Initial  Stress 

As  indicated  by  Bolotin  [7] , the  potential  per  unit  surface 

area  U associated  with  small  shell  motions  v,u,w  away  from  an 
go 

initial  equilibrium  state  with  stresses  o is  given  by  the 

jk 

quadratic  form 


U = 1 
go  2 


..h/2 

J -H 

2 


a / V « 1 

jk  \ j it 


V u \ dz 
. k i 


(12) 


where  the  V u are  the  displacement  gradients  in  Eqs.  (2),  and 

j 1 

summation  is  implied  on  the  subscripts  i,j,k.  After  substitution 
of  Eqs.  (1)  into  (2)  and  recognizing  that  Love's  postulate 


requires  O ~ o 

s z dz 

thickness  gives 


a - 0 , integration  through  the  shell 
z z 
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N 2 2 2 2 2 

U = _1  ■ V + V,.  + (sin  0 U + COS0  w)  +u,fl  + w , n 


2 2 
' r 


+ 2 v , ^ ( sin  0 u+cos  0 w) -2  v ( sin  0 u +cos0w,^)  J + 

oj  2 2 2 2 

+M  cos  0 v +v,0  - u,0  + w,0  - (sin  0 u + cos0  w) 


-2w(w,  + r s in  0 w,  ) - __2_  (v,  + sin0  u )(w,  + r sin0w,  ) 

00  s cos0[  0 *00  s 

+U,l9(W’s 0 r-sin<^  w»n)  + v(r  sin<£  w»;_^  " w’^)l  I 


s e e 


o 

2 2 2 

M 

N 

v , + u , + w ,1 

+ s 

2 v , 

s 

S S |s 

r 

* s 

- sin0(cos 0 v-w,„ ) -2rw,  u, 

u s s 

r s 


+ s 0 2v,  ( sin0u  + cos0w  )-  2v(  sin  0 u,  + cos0w,  ) 

r s 1 S s 


+2v,  v + 2w , w,  + u,  u,  + a u,  u, 
s 0 s 6 s 6 s 0 


(v  +)  h 


+ sin  0 u + cos  0 w)  ( -w,  + sin0  (w, 

0 Os  0 


COS0V)  + V,  (cos  0 V,  -w,  -r  sin  0 w,  ) 

j / s dee  s 

+ r w,  (sin  0 v - u,  ) - u,  (rw,  + sin  0 w,  ) 

ss  0 s s 0 0 

o 

I _ 1 M u,u,  cos  0 
2 s 0 s 0 
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Note  that  in  Eq . (13)  all  N and  M have  superscripts  (°) 

ij  ij 

which  denotes  the  initial  stresses.  In  utilizing  this  expression 

o e o e 

N , M = (14) 

ij  ij  ij  y 

e e 

where  N , M denote  the  stress  resultants  obtained  from  an 

ij  ij 

analysis  in  which  initial  stress  effects  are  ignored  and  is  a 
scale  factor . 

/ * 

II. 5 External  Potential  Due  to  Pressure  Loadings. 

As  discussed  by  Cohen  [8]  normal  pressure  loadings,  which 
are  follower  type  forces,  will  admit  to  a quadratic  potential 
only  for  certain  special  conditions  on  the  boundaries  of  the  shell. 
For  this  restricted  case,  the  external  potential  per  unit  undeformed 
surface  is  a functional  given  by  the  expression 


U =/9  • p_ 
eo  r 


1 1 *■ 

(v,  + s in  0 u+ru,  ) • w + — cos  0 (v  +w  ) 

Q s 2 


(15) 


e e 

The  pressure  p is  the  loading  which  is  used  to  generate  the  N 

e ij 

and  M in  Section  II. 4,  and,  as  shown  in  Figure  1,  is  assumed 
ij  e 

to  act  in  a direction  opposite  to  that  of  w.  In  this  form  p 

is  assumed  to  be  independent  of  s,0. 

II. 6 Hamiltonian  for  Shell. 

The  Hamiltonian  used  to  characterize  the  behavior  of  the 
conical  frustum  of  the  shell  is  the  integral  of  the  potentials 
in  Sections  II. 3-5  over  the  surface  of  the  shell, 

27T  _ s 


H = T-V  = 


If 

J n -/-s 


T - ( U + U 
o 

o go 


+ U + V ) 
eo  e 


rd  6 ds 


(16) 
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This  functional  gives  all  the  information  required  to  model  the 
shell.  Although  a set  of  differential  equations  and  associated 
boundary  conditions  could  be  derived,  in  order  to  develop  an 
algebraic  representation  of  the  structure,  a Fourier  series  is 
used  to  specify  the  dependence  of  the  shell  variables  in  the 
circumferential  angle  0 , and  an  explicit  polynomial  represen- 
tation is  used  to  specify  the  dependence  in  s. 

II. 7 Fourier  Analysis. 

Since  the  shell  is  closed,  i.e.,  0 < 0 < 2ir  all  variables 
required  to  describe  the  shell  may  be  expanded  into  a Fourier 
series  in  the  circumferential  angle  0 . In  particular 

N n N 

v (s , 0)  = v*°  (s)  - 52  v (s)sin  n0  + 52  v*ncos  n0 

n=l  n=l 

o N n N 

u(s,0)  = u (s)  + 52  u (s)cos  n6  + 52  u*nsin  n0  (17) 

n= 1 n= 1 

N n N 

w(s,0)  = w°(s)  + 52  w (s) cos  n 0 + ^ w*nsin  n 0 

n=l  n=l 

where  n is  the  circumferential  mode  number  and  N is  the  maximum 

number  of  Fourier  forms  considered  in  the  analysis.  Substitution 

of  Eqs.  (17)  into  Eqs.  (16)  and  integration  will  remove  the  angle 

0 from  the  functional,  at  the  expense  of  introducing  the  scalar 

parameter  n=l, ,N.  Fortunately,  the  portions  of  the  Hamiltonian 

associated  with  U and  T are  quadratic  in  all  the  field  variables 

o o 

so  that,  due  to  the  orthogonality  of  the  Fourier  series,  a complete 
decoupling  occurs  between  variable  sets  with  differing  mode  number 
n.  Further,  all  variables  with  asterisk  (*)  decouple  with  the 
variable  not  so  labeled,  so  that  a major  simplification  exists 


in  using  a Fourier  series  so  far  as  U and  T are  concerned. 

o o 

Unfortunately,  the  simplifying  features  associated  with 

the  quadratic  forms  for  U and  T do  not  carry  over  to  the 

o o 

more  complex  potentials  associated  with  initial  stresses,  and 

external  pressure,  since  in  these  potentials  three  field 

variables  interact,  two  displacement  gradients  and  an  initial 

stress.  Thus  axisymmetric  initial  stresses  interact  with  sets 

of  functions  for  all  mode  numbers  and  non-axisy'mmet  ric  initial 

stresses  will  couple  all  the  harmonics  together. 

In  order  to  attempt  to  develop  a useful  analysis  and  obtain 

a tractable  computer  program  only  the  initial  stresses  associated 

with  n=0  and  symmetric  (v=0)  motion  are  considered.  That  is, 

o o o o o 

only  the  initial  effects  associated  with  p , N , N , M and  M 

' e 6 s 0 s 

are  included  in  the  geometric  and  (quadratic)  external  potentials. 

e e e 

Torsional  prestresses  N , M and  M will,  for  a given  wave 

s Q s$n  s 6 

number  n,  couple  the  motions  u with  un*  and  so  will  not  be 

i i 

considered  at  this  time. 

Given  the  special  form  of  the  initial  stress  and  external 
potential  described  in  Section  II. 6,  then  a complete  decoupling 
of  the  shell  energies  is  possible  and  the  analysis  of  the  shell 
may  be  carried  out  one  mode  at  a time.  For  each  mode  the 
potentials  may  be  written  as  matrix  forms 


"n  ■^(|6JI[\]l5nl  + l4!|Ipn]l6n|)tdS 

T„ " 7 /(ntK1IMhCitpjiCi)-, 

+ Krp;]i5*i) ras 


u + u = 

gn  e n 


P 7T 
2 


(18a) 


(18b) 


(18c) 
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where  the  or' s are  time  dependent  functions. 


For  the  second.  Type  II  element 


Thus  the  shell  surface  displacements  and  displacement  derivatives 

5 [ may  be  expressed  through  the  meridionally  dependent  transfor- 
* n ' 

nation  A(s‘)  in  terms  of  the  polynomial  coefficients  |a  , defined 

. ' n ’ 

for  the  frustum  as  a whole 


1 r I1  i 1 

j 5 (s',t)|  = ^A(s’)  j j a ( t ) I 


A similar  expression  holds  for  the  Type  II  element.  The  time 
dependent  constants  may  in  turn  be  defined  in  terms  of  shell 
displacements  at  specified  positions  along  the  meridion  using 
Hermite's  Formula  [9].  For  the  Type  I element  the  locations 
chosen  are  s = -s  and  s = + s . At  each  location  i four 


quantities  are  defined,  u , v , w and  w . These  shell  dis- 

i i i i , s 

placements  may  in  turn  be  transformed  into  displacements 
n n _n  _ n 

u , v , w and /3  , termed  nodal  displacements  (Fig.  5)  which 

i i i i 

are  shell  displacements  in  the  axial,  circumferential,  radial 
(not  normal)  directions  and  rotation,  such  that 


n n n 

u = ”u  cos  <p  + w sin<£> 
i i i 

n n 

v = v 
i i 


w = -u  sin<£  + w cos  <f) 


w = (i 
i , s i 


The  results,  in  matrix  form  are  as  follows  for  the  Type  I element. 
In  this  matrix  C is  cos<f)  and  S is  sin0. 


In  Eqs.  (25,26)  nodes  i,j,k  are  at  s 


o 


= - s , 0 , and  s 
o 

respec  tively . 

The  finite  element  model  is  then  described  by  the  following 
energy  forms,  for  the  Type  I element: 


(27) 


Energy  expressions  for  the  Type  II  element  are  identical  except 
that  superscript  II  appears  in  place  of  I.  Equations  (27)  may 
be  simplified  to  read 


1 

1 
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(28) 
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gn  en  2 
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n 

matrices 

I ‘ 
K 

9 

■ r 

M 

and 

I ‘ 
G 

n _ 

. n 

n 

Similarly,  expressions  hold  for  Type  II  elements  in  Eqs.  (28) 

are  the  stiffness,  mass  and  geometric 
stiffness  for  the  shell  element.  Application  of  Hamilton's 
principle  gives  the  equations  which  characterize  the  vibration 
and  stability  characteristics  of  the  frustum  using  the  Type  I 
e 1 erne  n t , 
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¥ 


This  equation  has  eight  independent  degrees  of  freedom.  Equations 
similar  to  Eqs.  (29)  hold  for  the  Type  II  element,  but  have  twelve 
degrees  of  freedom. 

III. 2 Numerical  Aspects  of  Finite  Element. 

The  matrices  characterizing  stiffness  mass  and  shell  stability 
of  the  type  in  Eqs.  (27)  are  very  complex  functions  of  the  ordinate 
s.  Therefore,  these  matrices  are  evaluated  numerically  using 
Gaussian  quadrature  formulas,  with  variable  point  integrations 
depending  upon  the  degree  of  accuracy  required.  In  detail,  the 


’ I 

T 

r i 

' i • 

matrix  expressions  in  Eq.  (27),  e.g.. 

A 

K 

A 

are  first 

» 

L nJ 

. 

multiplied  and  evaluated  at  the  required  Gaussian  points,  weighting 
factors  are  then  multiplied  to  the  matrix,  and  the  result  is 
stored.  Results  are  accumulated  in  storage  as  the  program  pro- 
gresses through  the  various  integration  points.  It  should  be 
noted  that  since  the  element  matrices  are  symmetric  only  diagonal 
and  upper  tii angular  portion  need  be  calculated,  and  also  that 
a formal  matrix  product  may  be  replaced  with  much  simpler  and 
more  efficient  direct  statements  of  the  matrix  product. 

III. 3 External  Potential 

In  order  to  implement  the  conical  frustam  element  in  the 
solution  of  a finite  element  forced  response  problem  it  is  neces- 
sary to  develop  the  linear  external  potential  of  the  applied  loads. 
These  loads  are  either  statical  or  dynamical  and  are  distributed 
upon  the  surface  of  the  elements  and  around  the  nodal  rings  of 
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the  assemblage  ot  elements.  A stability  problem  is  riot  a linear 
forced  response  problem,  but  in  order  to  calculate  the  prestress 


state  it  is  first  necessary  to  perform  a static  analysis. 

Two  sets  of  applied  loads  are  defined.  The  first  set  is 

the  surface  loads  f , f and  f , assumed  to  be  applied  t.o  the 

s 0 £ 

middle  surface  of  the  shell,  permit  undeformed  surface  area,  in 
the  s,  6 and  £ directions,  respectively,  see  Figure  5.  (The  con- 
tribution to  the  external  potential  due  to  the  pressure  acting 
on  the  deformed  surface  area  is  accounted  for  in  Eq.  (15).) 

These  surface  loads  are  functions  of  s and  0 and  may  be 


represented  by 

the  Fourier 

series 

0 

N 

*n 

f (s 

,0)  = 

f (s) 

+ 2 

f 

(s)  cos  n0 

+ f sin  n0 

s 

s 

n=  1 

s 

s 

r 

o 

1 

N 

n 

*n 

f (s 

,0)  = 

f (s) 

-p  (s) 

+ 

2 

f (s) 

cos  n 6 + f 

sin  n0 

C 

5 

n=  1 

C 

C 

(30) 


*o  N 

f (8,0)  = f (s)  + I 

0 6 n = 1 


-f  (s)  sin  n0  + f cos  n0 
0 0 


The  second  set  of  applied  loads  are  line  loads  F , F , M and 

z r 0 

F are  simply  forces  and  a moment  per  unit  circumferential  length 

0 

which  are  applied  around  the  nodal  rings  and  are  represented  in 

terms  of  a Fourier  series  (F  , F and  M are  similar) 

z r (3 

o N n * n 

F (0)  = F + Z (F  cos  n0  + F sin  n0) 
z z n=l  z z 


'n 


o N n *i 

F (0)  = F + £ (F  cos  nO  + F sin  nO) 

r r n= 1 r r 

*o  N n *n 

= F + £ (-F  sin  n0  + F 

0 n = l 0 


0 


F (0) 

0 
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cos  n0 ) 


Eq.  (31)  continued 


M(0)=M  + £ ( M cos  n9  + H sinnd) 

3 3 n=l  3 3 


(31) 


The  external  potential  of  the  ring  loads  at  node  i denoted 

V is  then 
6i 

y'*2ir 

(F  ( 0 ) u (0)  + F w + Mft  + F v )r  d0  (32) 
o Zi  i ri  i Pi  °i  i 

Due  to  the  orthogonality  of  the  functions  in  0 the  forces  decouple 

> * 

into  harmonics  in  wave  number  n and  similar  expressions  for  the 
* th 

n series  at  node  i.  The  equivalent  nodal  load  for  the  n 
n 

harmonic,  (F  },  is  thus  given  as 
i 


-TTr 


r n— n 

F u + 


1 L 


Z 


n_m  n_n 

F w + M 3 + 

r 3 


n_n"j 
F v 
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Ji 


(33) 


The  external  potential  of  the  surface  loads  on  element  m is 

written  as  V is  expressed  by  the  integral 
m 


(s  , 9)  u (s  , 0) 


+ 


f w + 

C 


f V 

0 J 


r ( s ) ds 


d9 


(34) 


T 

Expressing  the  surface  loads  as  the  vector  {f}  = (f  f f ) and 

s C,  0 

using  the  element  transformations  defined  previously  the  surface 

n 

loads  may  be  converted  into  equivalent  nodal  loads  {F  } acting 

n m 

on  the  nodes  connected  to  element  m,  (q  },  as 

m 

= -tts  (q  } T f ^A(s’)J  {f  (s  ’ ) }r  (s ' ) ds1  (35) 

m o m ■'  - 1 

n T n 

= - (q  } {F  } 


As  implemented  in  the  computer  program  at  present  the 
integration  in  s'  to  compute  surface  loads  is  carried  out  in 
closed  form  for  surface  loads  assumed  to  be  constant  over  each 
element . 

III. 4 Assemblage  of  Global  Model 

Using  the  conical  frustum  element  it  is  possible  to  model 
any  shell  of  revolution  to  any  degree  of  accuracy  required  using 
a sufficient  number  of  elements.  Shells  having  meridional  curva- 
ture such  as  spherical  shells  will  require  relatively  more  ele- 
ments since  the  element  itself  lacks  meridional  curvature.  How- 
ever shells  composed  of  cylindrical,  conical  and  circular  disk 
components  may  be  modelled  with  very  few  elements,  especially 
when  using  the  higher  order  Type  II  element.  Assuming  the  global 

model  to  have  I nodes  then  the  system  will  have  four  degrees  of 

_n  _n  _n  _n 

freedom  per  node  (u  , w , 3 and  v ) or  41  degrees  of  freedom. 

i i i i _n 

The  system  degrees  of  freedom  are  assembled  in  the  vector  {q  }. 


System  stiffness,  geometric  stiffness  and  mass  matrices 


[c"j  and  [„"] 


are  formed  by  adding  up  the  element  matrices  taking 


into  account  the  nodal  connectivities  of  the  individual  elements. 

_n 

Similarly,  the  system  force  vector  {F  } is  formed  by  adding  up 
all  the  equivalent  nodal  loads  due  to  ring  loads  and  elemental 
surface  loads.  If  the  elements  are  assembled  sequentially  and 
there  are  no  branches  or  closed  loops  in  the  structure.  The 
resultant  equations  will  be  highly  banded.  Using  only  Type  I 
elements  the  bandwidth  will  be  only  8 degrees  of  freedom  and 
using  Type  II  elements  will  result  in  a bandwidth  of  12. 
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i (j»  t 

Assuming  a time  dependence  for  dynamic  motion  of  e the 

following  equations  of  motion  characterize  the  system. 


Using  this  equation  a variety  of  analyses  may  be  formulated. 

If  3 and  u)  are  set  equal  zero  the  static  analysis  is  obtained. 

From  the  initial  stress  state  thus  obtained  the  geometric  stiffness 
may  be  calculated.  Then  by  setting  F and  uj  to  zero  the  eigenprob- 
lem  for  the  buckling  parameter  3 is  obtained.  Setting  F to  zero 
and  specifying  3,  the  eigenproblem  for  the  vibrational  frequencies 
co  and  mode  shapes  is  obtained,  taking  into  account  the  destabil- 
izing effects  of  the  prestress.  If  3 and  u>  are  specified  the 
forced  vibration  problem  with  F a function  of  frequency  may  be 
solved.  Results  of  sample  cases  performed  are  presented  in  the 
next  section. 

IV  Implementation  and  Results 

In  order  to  obtain  results  from  the  methodology  presented 
here  a computer  program,  FASOR  (Fourier  Analysis  of  Shells  of 
Revolution),  has  been  developed.  The  program  has  been  designed 
to  minimize  core  requirements  and  computer  solution  time. 

Through  the  use  of  variable  dimensioning  of  all  subroutines 
the  amount  of  core  used  is  exactly  the  amount  needed  for  a par- 
ticular problem.  Sparse  matrix  techniques  have  been  employed 
in  most  matrix  manipulations.  Banded  matrix  techniques  have 
been  employed  in  the  static  analysis  back  substitutions  and  in 
the  sub-space  iteration  eigensolut ions . The  sub-space  iteration 
technique  [lO,ll]  is  very  efficient  and  calculates  only  the 
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lowest  buckling  and  vibrational  modes  desired  by  the  user. 

(Through  the  use  of  frequency  shifts  it  is  possible  to  solve  also 
for  either  the  rigid  body  modes  of  the  structure,  or  frequencies 
in  the  vicinity  of  a specified  value.) 

All  analysis  paths  oi  FASOR  have  been  checked  out  and 
results  compared  to  available  closed  form  and  numerical  solu- 
tions. Since  the  shell  theory  employed  herein  is  ,a  third  order 
theory  it  pertains  to  thicker  shells,  i.e.  larger  thickness  to 
radius  ratios,  than  do  most  available  codes.  Unfortunately, 
scarcity  of  available  solutions  for  thick  shells  has  so  far 
inhibited  the  investigation  of  FASOR's  solution  accuracy  compared 
to  other  programs.  This  and  the  investigation  of  the  effect  of 
shell  bending  stresses  on  stability,  which  is  also  unique  to 
FASOR,  will  be  pursued  in  future  studies.  In  order  to  make 
meaningful  comparisons,  the  cases  presented  in  the  following  are 
for  relatively  thin  shells  (h/a  ~ 0.01)  where  correspondence  with 
available  solutions  will  be  closest. 

Many  of  the  static  test  cases  may  be  obtained  from  Timoshenko 

[l2].  The  first  example  is  the  problem  of  a semi-infinite 

circular  cylinder  with  applied  end  shear,  H , and  moment,  M . The 

o o 

damped  sin-wave  solution's  results  for  end  radial  displacement 
and  rotation  for  a shell  having  the  following  properties: 

a = 10",  h = 0.1",  E = 30,000  Kips/sq.  in.,  V = 0.3;  are 

w ( 0 ) = .08569  H - .11015  M (37) 

o o 

B ( 0 ) = -.11015  H + . 28318  M 

o o 
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Usi.ig  ')  T/pe  II  clemi.nt'i  FA50R  gives 


3(0) 


o o 

1 100 A H + .28303  M 


(38) 


Using  only  one  Type  II  element  the  results  are  between  8%  and 
25%  in  error.  Various  circular  plate  bending  problems  have  been 
used  for  different  harmonics  and  loading  conditions.  For  these 
cases  displacement  results  were  exact  to  four  digits  or  more 
using  only  one  Type  II  element.  The  in-plane  stiffness  of  the 
elements  has  been  checked  out  using  various  closed  form  solutions 
for  conical  membrane  shells  under  loading  conditions  such  as 
internal  pressure,  symmetric  and  anti-symmetric  weight  loading  and 
ring-loads.  All  results  were  very  accurate  using  not  more  than 
five  Type  II  elements. 

A cylinder  vibration  problem  was  compared  to  Forsberg's 

results  [l3]  using  Fliigge's  cylindrical  shell  theory.  The  follow- 

6 

ing  dimensions  and  properties  were  used:  E = 30x10  psi, 

2 

v = 0.32,  p = .000253886  lb-sec  /in.,  a = 3.045  in.,  h = .06  in., 

length  “ 16.2  in.  For  classical  simply  supported  boundary  conditions 

w=v=N  = M =0in  wave  n = 4 Forsberg's  lowest  frequency  is 
s s 

974  cy/sec.  Using  10  Type  II  elements  FASOR  gave  973.47  cy/sec. 

A comparison  for  the  geometric  stiffness  effects  of  the 

finite  element  is  provided  by  the  buckling  problem  of  a conical 

frustum  loaded  by  external  pressure  [I*].  The  shell  investigated 

6 

has  the  following  properties:  E = 30x10  psi,  V = 0.3,  h - 0.1 
in.,  left  radius  = 5 in.,  right  radius  = 10  in.,  cone  angle 
0 = 30°,  length  2s  = 10  in.  The  boundary  conditions  are  at  each 

25 


, 


0. 


The  lowest  (critical)  buckling  pressure 


end , u = w = M = N 

s s 9 

was  found  for  circumferential  wave  n = 9 and  is 


Pcr  = 

80.53 

p s i 

: Re  f 

. 13 

L 

Pcr  = 

81.99 

p s i 

: 10 

Type 

II 

elements 

Pcr  = 

82.72 

p s i 

: 5 

Type 

11 

elements 

Pcr  = 

98. 

ps  i 

: 1 

Type 

II 

element . 
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